//
//  Move.cpp
//  phi_four
//
//  Created by Snir Gazit on 3/10/12.
//  Copyright (c) 2012 __MyCompanyName__. All rights reserved.
//

#include "Move.h"
// General Move

void Move::init(MTGen* eng, Lattice * lat) {
	eng_ = eng;
	lat_ = lat;
	zeroOne_.init(eng, UniRealParams(0, 1));
}
bool Move::accept(double p) {
	if (p > 1)
		return true;
	double draw = zeroOne_.draw();
	if (p > draw)
		return true;
	else
		return false;
}
bool Move::checkMove() {
	int totalBonds;
	return lat_->testConfiguration(totalBonds);
}


void CreateLoopMove::init(MTGen *eng, Lattice *lat,
		QmcParameters &qmcParameters, BoltzWeightType boltzRatio) {
	static_cast<Move *>(this)->init(eng, lat);
	chooseSite_.init(eng, UniIntParams(0, lat->getNumOfSites() - 1));
	setCoupeling(boltzRatio,qmcParameters.g.getVal());
}
void CreateLoopMove::setCoupeling(BoltzWeightType boltzRatio,double g) {
	g_ = g * 0.5;
	RemoveRatio_[0] = 0.;
	for (int j = 0; j < MAX_BOND; j++) {
		CreateRatio_[j] = boltzRatio[j] * g_;
		if (j != (MAX_BOND - 1))
			RemoveRatio_[j + 1] = 1 / (boltzRatio[j] * g_);
	}
}
// Create Loop Move

int CreateLoopMove::makeMove() {
	int wormLength = 0;
	currentLoop_.clear();
//currentLoop_.resize(avWormLength_);
// Choose First Site
	int jumpSite = chooseSite_.draw();
	Lattice::SiteIter jump = lat_->getSite(jumpSite);
	lat_->setIra(jump);
	lat_->setMasha(jump);
	getWeights();
	double normFirst = accumulate(weights_.begin(),weights_.end(),0.);
// Start Looping
	do {
		wormLength++;
		nextStep();
		getWeights();
	}while (!lat_->isClosed());
	double normLast = accumulate(weights_.begin(),weights_.end(),0.);
	double pc = normFirst / normLast;
	if (!accept(pc)) {
		eraseLoop();
		eraseWorm_(1);
		//	assert(this->checkMove());
		return 0;
	}
//	assert(this->checkMove());
	eraseWorm_(0);
	return wormLength;
}
double addRatio(int totalBonds, int otherTotalBonds) {
	return double(totalBonds + 1) / double(2 + totalBonds + otherTotalBonds);
}
double removeRatio(int totalBonds, int otherTotalBonds) {
	return double(1 + totalBonds + otherTotalBonds) / double(totalBonds);
}
void CreateLoopMove::eraseLoop() {
	Lattice::SiteIter cIra = lat_->getIra();
	for (auto itl = currentLoop_.begin(); itl != currentLoop_.end(); itl++) {
		auto &cBondData = *(itl->bondData);
		Lattice::SiteIter nextIra = cBondData.neighbourSite_;
		// int &totalBonds = nextIra->totalBonds_;
		int &bondVal = cBondData.getBondVal(loopType_);
		if (cBondData.createRemove_ == CREATE) {
			bondVal--;
			(nextIra->getTotalBonds(loopType_))--;
		} else {
			bondVal++;
			(cIra->getTotalBonds(loopType_))++;
		}
		cIra = nextIra;
	}
}

void CreateLoopMove::getWeights() {
	auto itW = weights_.begin();
// Choose First Site
	auto ira = lat_->getIra();
	auto &cBonds = ira->XYBondArray_[loopType_];
	int cTotalBonds = (ira->getTotalBonds(loopType_)), cOtherTotalBonds =
	ira->getTotalBondsOther(loopType_);
	//currentLoopData_.totalBonds = cTotalBonds;
	double weightRemove;
	if (cTotalBonds == 0)
	weightRemove = 0;
	else
	weightRemove = (RemoveRatio_[cTotalBonds + cOtherTotalBonds]
			* removeRatio(cTotalBonds, cOtherTotalBonds));
	for (auto itb = cBonds.begin(); itb != cBonds.end(); itb++) {
		auto bond = *itb;
		int bondVal = bond.getBondVal(loopType_);
		if (bond.createRemove_ == CREATE) {
			Lattice::SiteIter nSite = bond.neighbourSite_;
			int totalBonds = nSite->getTotalBonds(loopType_);
			int otherTotalBonds = nSite->getTotalBondsOther(loopType_);
			double weight = (CreateRatio_[totalBonds + otherTotalBonds]
					* addRatio(totalBonds, otherTotalBonds))
			/ double(bondVal + 1);
			*itW++ = min(1., weight);
		} else {
			double weight = weightRemove * double(bondVal);
			*itW++ = min(1., weight);
		}
	}
}
int CreateLoopMove::chooseDirection() {
	double totalWeight = accumulate(weights_.begin(), weights_.end(), 0.);
	//currentLoopData_.totalWeight = totalWeight;
	double draw = totalWeight * zeroOne_.draw(), sum = 0.;
	int chosen = -1;
	for (auto itW = weights_.begin(); sum < draw; itW++) {
		++chosen;
		sum += *itW;
	};
	if (chosen == -1)
	chosen = 0;
	assert(!(chosen < 0));
	assert(!(chosen >= weights_.size()));
	return chosen;
}
void CreateLoopMove::nextStep() {
//	chooseDirection_.setParams(DiscreteVariableParams(weights_));
	int chosenDir = chooseDirection();
	Lattice::SiteIter cIra = lat_->getIra();
//    int spaceDir=chosenDir/(2*DIM);
	auto &cBondData = cIra->XYBondArray_[loopType_][chosenDir];
	Lattice::SiteIter nextIra = cBondData.neighbourSite_;
// int &totalBonds = nextIra->totalBonds_;
	int &bondVal = cBondData.getBondVal(loopType_);
	if (cBondData.createRemove_ == CREATE) {
		bondVal++;
		assert(bondVal < MAX_BOND);
		(nextIra->getTotalBonds(loopType_))++;
	} else {
		bondVal--;
		assert(bondVal >= 0);
		(cIra->getTotalBonds(loopType_))--;
	}
	lat_->setIra(nextIra);
	currentLoopData_.bondData = &cBondData;
	currentLoop_.push_back(currentLoopData_);
//	assert(this->checkMove());
}

double CreateLoopMove::getEraseProb() {
	return boost::accumulators::mean(eraseWorm_);
}

void CreateLoopMove::setAvWormLength(int avWormLength) {
	avWormLength_ = avWormLength;
}
int CreateLoopMove::getLoopType() const {
	return loopType_;
}
void CreateLoopMove::setLoopType(int loopType) {
	loopType_ = loopType;
}
void CreateLoopMove::switchLoopType() {
	loopType_ = (loopType_ + 1) % 2;
}
